package scu.maqiang.piezoelectricity;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L2PE;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class APE1DL2 {

	public static void main(String[] args) {
		Mesh1L2 mesh = new Mesh1L2();
		mesh.line(0, 1, 100);
		mesh.translate(1.0);
		
		FES1L2PE fs = new FES1L2PE(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		
//		double G1 = 1.0, ep1 = 0.5, e1 = 0.1;
//		fs.assemblePEStiff(new double[] {G1, ep1, e1}, BVPType., A);
		
		double G1 = 1.0, G2 = 2.0, lambda = 0.5, ep = 0.2, e1 = 0.4, e2 = 0.3;
		fs.assemblePEStiff(new double[] {G1, G2, lambda, ep, e1, e2}, BVPType.AXISYMMETRIC, A);
		double f1 = 1.0, rhof = 2.0;
		fs.assemblePESource(new double[] {f1, rhof}, BVPType.AXISYMMETRIC, RHS);
		double[] x = new double[fs.GetNdof()];
		fs.applyBC_MBN(A, Direct.All, 1);
		fs.applyBC_MBN(RHS, Direct.All, 0.0, 1);
		

		
		NewIterSSolver solver = new NewIterSSolver(A);
		//solver.GMRES_SSOR(RHS, x, 1.5, 20, 1);
		//solver.PCGSSSOR(RHS, x, 1.5, 2);
		Arrays.fill(x, 0.0);
		//solver.CGNE(RHS, x, 1);
		//solver.PCGNE_M(RHS, x, 1);
		System.out.println("x min:" + MVO.min(x) + "max: " + MVO.max(x));
		Arrays.fill(x, 0.0);
		solver.CGNR(RHS, x, 1);
		System.out.println("x min:" + MVO.min(x) + "max: " + MVO.max(x));
		
		Arrays.fill(x, 0.0);
		DirectSSolver directSolver = new DirectSSolver(A);
		x = directSolver.LUSolve(RHS, 1, 1e-6);
		System.out.println("x min:" + MVO.min(x) + "max: " + MVO.max(x));

		double[][] u = new double[4][mesh.getNv()];
		fs.extractComponent(x, u);
		
		//ѹ�����⾫ȷ��
//		ScalarCoefFunc uExact = (xy, label, param) -> 
//		                        (rhof * e1 - f1 * ep1)/(G1 * ep1 + e1 * e1) * (0.5 * xy[0] * xy[0] - xy[0]);
//		ScalarCoefFunc phiExact = (xy, label, param) -> 
//		                        -(f1 * e1 + rhof * G1)/(G1 * ep1 + e1 * e1)* (0.5 * xy[0] * xy[0] - xy[0]);
		
//		u[2] = fs.valueFromFunc(uExact, null);
//		u[3] = fs.valueFromFunc(phiExact, null);
		mesh.toTecplot("APESolution.dat", u);		
	}
	
	
}
